Tetrahedral clustering in molten lithium under pressure 
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A series of electronic and structural transitions are predicted in molten lithium from first prin- 
ciples. A new phase with tetrahedral local order characteristic of sp^ bonded materials and poor 
electrical conductivity is found at pressures above 150 GPa and temperatures as high as 1000 K. 
Despite the lack of covalent bonding, weakly bound tetrahedral clusters with finite hfetimes are 
predicted to exist. The stabihzation of this phase in lithium involves a unique mechanism of strong 
electron localization in interstitial regions and interactions among core electrons. The calculations 
provide evidence for anomalous melting above 20 GPa, with a melting temperature decreasing below 
300 K, and point towards the existence of novel low-symmetry crystalline phases. 
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The properties of materials under high pressure can 
change in profound and unexpected ways. At ambient 
pressure (P) and temperature (T), Li can be regarded as 
the prototype for an ideal metal. However, theoretical 
studies [H have predicted that as a result of increased 
core-valence electron interactions, low-T dense solid Li 
would undergo a series of symmetry-breaking transitions, 
culminating in a Li2-paired crystal with semi-metallic 
properties at pressures above 150 GPa. Besides the regu- 
lar transition at 7.5 GPa from body-centered cubic (bcc) 
to face-centered cubic (fee) structure measurements 
to date 0] have confirmed only an initial transition at 
39 GPa from fee to a less compact structure with a 16- 
atom cubic cell (c/16 space group) via an intermediate 
disordered phase (hRl). The existence and properties 
of even lower-coordinated structures at higher pressures 
remains an open question due to the lack of experimen- 
tal data and difficulties associated with conclusive the- 
oretical predictions of low-temperature thermodynamic 
crystalline stability. 

It was shown recently 4] that solid-solid transitions in 
Na are preceded by analogous changes in its liquid but 
at much lower pressure. Similar behavior was also sug- 
gested for other light alkalis. This raises the intriguing 
question, answered in what follows, of whether the ex- 
istence of low-symmetry structures in Li d i, i, [1, 
could be established by investigating its molten phase. 
Molecular dynamics simulations of liquids have the ad- 
vantage of depending minimally on initial conditions and 
are not subject to the inherent bias of comparing a finite 
number of solid configurations. Finding low-coordinated 
structures in molecular dynamics simulations of a liquid, 
in addition to being remarkable in itself, would therefore 
constitute rather conclusive evidence for the existence of 
similar crystalline phases. 

We have carried out first principles molecular dynamics 
(FPMD) simulation of liquid Li in the density range cor- 
responding to 3.06 > Ts > 1.60 (where |7r(rsao)^ = V/N, 
ao is the Bohr radius, V the volume and N the number of 



valence electrons) at temperatures up to 3000 K, and of 
the known crystalline phases: bcc, fee, and c/16, between 
and 90 GPa. Density functional theory with a plane 
wave basis set and the Born-Oppenheimer approximation 
were used for all simulations |8j]- They were performed 
in the NVT ensemble with cubic supercells and periodic 
boundary conditions. For the simulations of bcc and c/16 
solids, and all liquids we used 128 atom supercells with F- 
point sampling of the Brillouin zone. Simulations of the 
fee solids were performed with 108 atoms and 8 k-points 
(the F-point only was found to be insufficient). Most 
simulations ran for at least 10 ps, with some up to 20 
ps, while for obtaining the melting temperatures we have 
computations as long as 200 ps. Extensive convergence 
tests were carried out for size effects, various simulation 
parameters, and the validity of the pseudopotential ap- 
proximation Q . Good agreement was also obtained with 
existing liquid experimental data [l3| . 

Results for the structural analysis of the liquid along 
a 1000 K isotherm are presented in Fig. [TJ Initially, the 
first peak of the pair correlation function, g{r), broad- 
ens (Fig. [Ija)) in a way previously observed in Na Q]. 
Upon further compression, it splits entirely, indicating 
significant further lowering in the coordination. For a 
more detailed analysis of the structure, we examine the 
evolution of neighbor distances with increasing density 
(Fig.[ljb)). While there is no indication for discontinuous 
liquid- liquid transitions, several distinct regions with dif- 
ferent liquid structures can be identified, which correlate 
well with the electronic and melting properties discussed 
further below. These are: (i) r, > 2.60 {P < 23 GPa); 
(ii) 2.60 > > 2.05 (23 GPa < P < 150 GPa); and 
(in) rs < 2.05 (P > 150 GPa). The initial changes are 
analogous to those observed in liquid Na [3] - a transition 
from a bcc- like to an fee-like local order in (i) , followed in 
(ii) by lowering in the coordination (number of neighbors 
under the symmetrized first peak of g{r)) to 8 + 4 + . . ., 
and the liquid acquiring a c/16-like local order. 

In (iii) , the density- rescaled average distances to neigh- 
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FIG. 1: Structural changes in molten Li along the 1000 K 
isotherm, (a) Pair correlation function, g(r), for selected den- 
sities. Vertical dashed lines indicate the average positions of 
the 4*^, S*'', and 12*'^ neighbors, (b) Evolution of density- 
rescaled average interatomic distances with pressure and den- 
sity. First, average interatomic distances are computed to the 
first (ri), second (r2), and so on, neighbors. Next, density- 
rescaled distances are calculated as di = {ri/rs)/{ri /r^), 
where — 3.06 and refers to distances obtained at r^. Fi- 
nally, dm-n is the averaged value of di for neighbors, i, from 
m to n. The corresponding values of m and n are indicated 
next to each curve. The inset shows bond angle distributions 
between the first two neighbors in molten Li (solid curves) for 
Vs = 3.06 (blue) and 1.80 (red), of solid bcc Li (dashed blue 
line) at 500 K and = 3.06, and of molten diamond [ll[ 
around 250 GPa (dashed red curve). Angles among nearest 
neighbors in ideal bcc (70.50°) and diamond (109.50°) crys- 
tals are indicated. The bcc data are scaled by a factor of 1/3 
relative to the other angle distributions. 



bors 1-4 continue to contract even faster, while the next 
eight, as well as 13-16, move away. As a result, by 
rs = 1.80 (P ~ 370 GPa), the first coordination shell 
completely splits in two, the coordination becomes only 
four and remains roughly so up to at least — 1.6 
(P ~ 810 GPa). The distribution of angles among the 
nearest neighbors also becomes rather unexpected (inset 
in Fig. [TJb)) as it has a peak at 109.5°. The parallel 
with the liquids of materials with sp'^ bonding is indeed 
striking when compared to the corresponding angle dis- 
tribution in molten carbon, obtained by melting diamond 
at similar pressure (inset in Fig. [Ijb)). Therefore, there 
is a large pressure range, 150 GPa < P < 810 GPa, 
for which we predict that Li has tetrahedral local order, 
hitherto not seen in a liquid metal, but characteristic of 
semiconductor liquids with sp'^ bonding. 

The apparent structural similarity with covalently 
bonded liquids has prompted us to investigate the persis- 
tence probability of Li clusters. For this purpose, we eval- 
uate a function Pn{t), which gives the statistical proba- 
bility for the first (n—l) nearest neighbours found around 
an atom at time t to be the same as those found at time 
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FIG. 2: Survival probability, Pn{t), of clusters as a function 
of their size, n. Results are shown for Li and for a hard sphere 
(HS) liquid. The times at which the P„(t)'s are evaluated are: 
43 fs (for rs = 3.06), 27 fs (r, = 2.22), 20 fs (r, = 2.00), 15 fs 
{rs = 1.80) and 12 fs (r^ = 1.60) for Li, and 20 fs for the HS 
liquid. The HS radius is chosen such that the self-diffusion 
constant and the short-distance (/(r)'s of the HS and Li liq- 
uids are roughly the same. Computations with HS's at lower 
densities show rapidly decreasing values of Pn{t); already an 
order of magnitude lower at = 2.0. The inset shows sur- 
vival probability of clusters as a function of temperature at 
rs = 1.80. 



t = 0. For completely uncorrelated particles, P„(i) would 
drop monotonically with increasing n for any fixed t. 
This is a result of the fact that the larger the cluster the 
more particles there are that can leave the cluster during 
any time interval. For the same reason, Pn{t) should be 
a decreasing function of n for a sufficiently large t even if 
there is some weak binding among the atoms. For exam- 
ple, if p is the average probability for a given particle to 
leave a cluster during a time interval t, then the survivor 
probability of that cluster is P„(t) ~ {\ — p){n — 1) for 
large t (and neglecting particles re-entering the cluster). 

In order to study cases with weak metastability, the 
time chosen for the evaluation of P„(i) should be suffi- 
ciently short such that results are non-trivial. The values 
oit at which we have evaluated Pn{t) for each density are 
listed in the caption of Fig. [51 The ratios among them are 
roughly equal to the ratios of the periods of atomic vi- 
brations at different densities, which have been estimated 
from a velocity autocorrelation analysis. This allows us 
to quantitatively compare the metastability of clusters 
for a range of densities, over which diffusion rates and 
vibrational properties differ significantly. As shown in 
Fig. [21 Pn{t) changes qualitatively when density is in- 
creased. The appearance of a peak at n = 5 is evi- 
dence for a weak metastability of the tetrahedral clusters. 
To estimate the importance of caging effects 1^, which 
could appear in a compressed liquid, we have also carried 
out simulations of hard-sphere liquids. The comparison 
with Li (Fig. [2]) demonstrates that the peak in P„(t) of 
Li is much higher, indicating significantly stronger inter- 
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FIG. 3: Electronic properties of dense molten Li. (a) Elec- 
tronic density of states, DOS. The corresponding ra val- 
ues for each density are listed in the legend, (b) Valence 
electron bandwidth (blue triangles) and integrated DOS of 
core electrons with p angular momentum character, p-DOS 
(red circles), (c) Calculated dc conductivity, a, along the 
1000 K and 3000 K isotherms. Error bars are statistical 
standard deviations. Comparison is made with experimen- 
tal measurements [l5| . The data from [ij are divided into 
two sets, representing two temperature intervals. Only quali- 
tative comparison is meaningful with [ij as a model equation 
of state has been used there to obtain Li density, T, and a. 



atomic correlations. 

We now examine the changes in the electronic prop- 
erties that are likely to drive the structural transitions. 
With increasing density, the DOS at the Fermi level grad- 
ually decreases ( Fig. [3]^a)). A similar effect, interpreted 
as a Peierls symmetry breaking, was observed in liquid 
Na except that its strength in Li is stronger, which 
can be understood in terms of an increasing hardness 
of the effective repulsive potential; it is similar, for ex- 
ample, to the decreased amplitude of the distortion in 
crystalline Sb in comparison to crystalline As 
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interpret these changes and their likely consequences, we 
first look at the valence electron bandwidth (Fig. [H^b)). 
ft increases with increasing density, es expeced upon den- 
sification. This tendency remains up to P ~ 23 GPa 
(ts = 2.6), which matches well the pressure range over 
which the liquid becomes more compact. Above this pres- 
sure, the bandwidth begins to decrease due to a develop- 
ment of partial (p-character) bonding, which lowers the 
band-structure energy. The DOS develops a marked peak 
well below the Fermi level. This is also the range over 
which the coordination in the liquid begins to decrease. 
Finally, starting at P 150 GPa (r-g < 2.05), the valence 
band begins to broaden again, but the DOS at the Fermi 
level does not decrease much further. This density range 
coincides with the conditions under which the tetrahedral 
coordination develops. 

The consideration of valence electrons alone is clearly 
insufficient to understand the modifications occurring at 
P > 150 GPa, which can be explained by examining the 



core electrons. Indeed, for < 2.05 (P > 150 GPa) 
the p angular momentum character of the core electrons 
increases rapidly (Fig. EJb)), indicating core-core over- 
lap. At the same time, the volume available to the va- 
lence electrons, now squeezed into interstitial regions, de- 
creases linearly with increasing density and hence the 
rapid increase of the valence bandwidth. The resulting 
"anti-diamond" structure (inset in Fig. [5]) is with ions 
forming a tetrahedral network, but the valence electrons 
occupying preferentially the voids of the diamond struc- 
ture instead of being located between adjacent ions. 

The electrical conductivity (u) of compressed liquid Li 
exhibits exactly the opposite P and T dependence usual 
for metals. Our computations [3] show that along the 
1000 K isotherm it has nearly a 10-fold drop between 3 
and 150 GPa. When heating the liquid to 3000 K near 
and above 250 GPa, there is about twofold increase in a, 
but this increase vanishes at lower P. These results are 
compared with measurements 14, 15| for which we now 
provide an explanation that is very different from what 
was originally proposed. 

First, the electron localization and decrease of DOS 
at the Fermi level lead to a significant drop in a with 
P. Second, the above-mentioned electronic changes are 
closely related to structural transitions from a higher to 
a lower-coordinated liquid. This trend is reversed when 
heating the liquid at a constant V; it reverts to a more 
homogeneous local order (see Fig. |4]) as favored by en- 
tropy. The result is an increase of a with T, as seen in 
our calculations above 150 GPa. This effect is of course 
countered by increased electron-ion scattering at high T. 
At P ~ 100 GPa, the two effects cancel, and at lower 
P where the structural changes are less significant, the 
scattering effects dominate. This explanation for the ob- 
served changes in a is different from what was previously 
suggested and we emphasize that the increase in conduc- 
tivity observed at the highest-P measurements by Fortov 
et al. [3]) is not a pressure but a temperature effect. De- 
spite the fact that the latter data were obtained using a 
model equation of state to estimate the density, resistiv- 
ity, and T, the general trend of increased a with T at 
sufficiently high P is consistent with our findings. 

Finally, we compute the melting curve (Fig. [5]) using 
a heat-until-melting approach, which provides an upper 
bound for the melting temperatures (Tm)- The method 
was shown [4j to be valid for Na for reasons that are likely 
to hold for Li. The melting curve has a steep negative 
slope between roughly 20 and 80 GPa, with drop- 
ping from (680 ± 25) K to (275 ± 25) K. The shape 
of the melting curve above bcc is relatively flat, which 
is consistent with a gradual transformation in the liq- 
uid from 6cc-like to /cc-like local order. The fee solid is 
denser than bee and the melting slope above it initially in- 
creases, in accordance with the Clapeyron equation. The 
onset of symmetry breaking transitions lead to lowering 
of the liquid electronic band structure energy, its densifi- 
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FIG. 4: Evolution of density-rescaled average interatomic dis- 
tances with pressure and density at T = 1000 K (left panel) 
and T = 3000 K (right panel). As expected, increasing tem- 
perature leads to more homogeneous structures. 



cation, and, hence, to the turnover of the melting curve 
above the fee phase. The maximum Tm is near 20 GPa 
- exactly where we identified the onset of the changes 
towards lower coordination. The anomalous melting be- 
havior persists until pressures where the solid also begins 
to undergo Peierls symmetry breaking transitions. 

We have only computed the melting curve in the pres- 
sure range where the crystalline structures are well estab- 
lished. Above 100 GPa, (or even down to 70 GPa), they 
remain unsettled; previous suggestions include Cmca [I] 
and Cmca24 We propose that future investigations 
of low-T crystalline phases of compressed Li focus on 
diamond-like tetrahedral structures. The appearance of 
such a phase, and especially its persistence at high T, 
is completely unexpected. While the modifications in Li 
are initially driven by Peierls-like distortions, at higher 
P they are determined by core-core electron interactions 
and valence electron localization (resulting from core- 
valence interactions). This behavior describes a distinct 
regime, likely present in other materials, where both va- 
lence and core electrons are responsible for chemical and 
physical properties. Such effects could have far reach- 
ing consequences in areas ranging from planetary mod- 
elin g to the study of superconductivity under pressure 



[lil bd. 21, [2^, [23|. Another interesting aspect is the 
possible implications of these results for the properties 
of dense hydrogen. It has been discussed |2J] that 
similarities between Li and H could be used to predict 
high-pressure phases, including superconductivity j25j, in 
metallic H. 
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FIG. 5: Mel ting curve of Li under pressure. Available experi- 
mental data axe shown in green. Uptriangle and down- 
triangle pairs indicate FPMD simulations of solid and liq- 
uids phases, respectively, on the same isochore, which bracket 
the melting temperature. The dashed line connecting the 
highest-pressure experimental data to the theoretical points 
near 10 GPa is only a guide to the eye. The vertical dotted 
lines indicate ambient temperature phase boundaries between 
solid phases. Only the high temperature phase boundary be- 
tween fee and c/16 has been estimated based on their re- 
spective melting temperatures. For pressures above 90 GPa, 
previous theoretical studies have proposed crystals with the 
Cmca 1] and Cmca24 0] symmetries; hence, we have indi- 
cated this phase as CmcaX. The inset shows the tetrahedral 
structure, which develops in the liquid above 150 GPa. The 
picture is from a liquid configuration taken at ~ 1.8 and 
T — 1000 K. The ions are shown as red balls, each connected 
to its nearest neighbours. Isosurfaces of the valence charge 
density are shown in gold; they illustrate the localization of 
valence electrons in the interstitial regions. 
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